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Abstract 

Computational mathematics plays an increasingly important role in 
computational fluid dynamics (CFD). The aeronautics and aerospace re¬ 
search community is working on next generation of CFD capacity that is 
accurate, automatic, and fast. A key component of the next generation of 
CFD is a greatly enhanced capacity for mesh generation and adaptivity of 
the mesh according to solution and geometry. In this paper, we propose 
a new method that generates triangular meshes on domains of curved 
boundary. The method deforms a Cartesian mesh that covers the domain 
to generate a mesh with prescribed boundary nodes. The deformation 
fields are generated by a system of divergence and curl equations which 
are solved effectively by the least square finite element method. 

1 Introduction 

Despite considerable success, mesh generation and adaptation remain a challeng¬ 
ing task for numerical simulation on complex geometries. A NASA sponsored 
study entitled CFD Vision 2030 states that Mesh generation and adaptivity 
continue to be significant bottlenecks in the CFD workflow, and very little gov¬ 
ernment investment has been targeted in these areas. As more capable HPC 
hardware enables higher resolution simulations, fast, reliable mesh generation 
and adaptivity will become more problematic. Additionally, adaptive mesh 
techniques offer great potential, but have not seen widespread use due to issues 
related to software complexity, inadequate error estimation capabilities, and 
complex geometries. 

In this paper, we propose an innovative approach to the generation and 
adaptivity of triangular meshes that overcomes problems of current techniques. 

We demonstrate the proposed method on a domain Z), which resembles 
the right half of the curved domain in Figures 8.1 of the well-known textbook 
[George 1991] Tj. In Figures 8.2, the so-called superposition-deformation method 
is described and demonstrated. This method, investigated by [Yerri et al. 1984] 
[2], [Cheng et al. 1988] 0 and [Shephard et al. 1988] [2], constructs a mesh 

1 Department of Mathematics, University of Texas at Arlington, Arlington, TX, US 

2 Nanhua Power Institute, Zhuzhou, China 

3 Department of Mathematics, Dallas Baptist University, Dallas, TX, US 


1 



of a curved domain with prescribed boundary nodes. A modified version is 
described in detail which uses quadtree-based partitioning of the squares near 
the boundary to better approximate the geometry of the boundary. After the 
squares near the boundary are sufficiently refined, all squares that intersect 
the domain are split into triangles. The boundary squares are treated based 
on a careful classification of the boundary patterns. The interior triangles are 
smoothed by iterations based on barycentric coordinates of the vertices. The fi¬ 
nal mesh generated by the modified method is shown in Figure 8.8. It is pointed 
out in [George 1991] [1] that (1) Treatment of squares near the boundary is a 
difficulty of the method; (2) The generated mesh respects the initial contour; 
but, in general, it may contain a slightly different boundary discretization in 
that a given edge being the union of smaller edges. In fact, visual inspection 
of the mesh in Figure 8.8 reveals poor quality of the mesh near the boundary. 
In contrast, our method deforms a Cartesian mesh that covers the domain to 
generate a mesh on the domain using only the prescribed boundary nodes. The 
deformation fields are generated by a system of divergence and curl equations 
which are solved effectively by the least squares finite element method. The 
deformation field preserves the Jacobian determinant and thus the method pre¬ 
vents element inversion. Moreover, the implementation of our method is based 
on differential equations, which do not rely on treatment of boundary squares 
on case-by-case bases. The resulted mesh uses the initial edges on the boundary. 
Thus the initial description of the geometry is preserved. In the next section, 
we use an example to demonstrate our method. In the following sections, the 
mathematical algorithm of the deformation method and the least squares finite 
element method used in the method are further discussed. 


2 An Example 


A uniform Cartesian mesh on the background domain [1, 9] by [1,12] is shown in 
Figure[la]below. The domain D is the interior of the contour described by points 
#1 to #18 on the curved boundary, and by the nodes (1, 5), (1,6), (1, 7), (1, 8), 


(1,9), and (1,10) on the vertical boundary, see Figure la We select a set of 18 


nodes on the Cartesian mesh that are close to the boundary of domain D 1 and 
move them to the points #1 through #18, respectively (see Appendix). For 
instance, by our method, node Pi(l,4) is moved to #1; P2(l,3) to #2; P 3 (2,2) 
to point #3; P 4 (3,2) to #4; node P 5 (4,2) to #5; P 6 (4,3) to #6; P~ 7 ( 5,3) to 
#7; Pg(6,4) to #8; node Pg(7,4) to point #9; node Pio(8, 5) to #10; Pn(7, 6) 
to #11; Pi 2 (6,7) to #12; P 3 (5,8) to #13; Pi 4 (4,8) to #14; Pi S (4,9) to #15; 
Pi6(3,9) to #16; Pn(2, 10) to #17; Pis(l, 11) to #18, etc. Correct movements 
of all other nodes are enabled by solving a divergence-curl system with suitable 
Dirichlet conditions on these 18 nodes as if they were boundary points. The 
selected background nodes are the small dots in Figure b, which are moved 
to the prescribed boundary nodes #1 — #18, respectively. In Figure [2a] the 
prescribed boundary of domain D is colored in blue. Finally, a triangular mesh 
is generated by connecting suitable pairs of the nodes on domain D , see Figure 
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Figure 2 


3 The Deformation Method 

We now describe the original deformation method in differential geometry, which 
we adopted to generate non-folding meshes. 
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3.1 The Deformation Method of Jurgen Moser in Differ¬ 
ential Geometry 

The method we propose has its origin in differential geometry. It was used in 
[Moser 1965] [5] for a study of volume elements in Riemannian manifolds and 
was modified for domains in R d (d is the dimension of the space) with boundary 
in [Dacorogna et al. 1990] [6]. They showed how to construct a diffeomorphism 
(i.e. a smooth and invertible deformation) from a 2- or 3-D domain 12 1 to D 2 with 
prescribed Jacobian determinant. The boundary nodes of 12 1 are matched to 
suitable boundary nodes of D 2 in a one-to-one correspondence. For a prescribed 
positive size function f(x, y , z) > 0, a method based on Poisson equations is 
devised which determines a suitable velocity vector field V(x,y,z,t). Each in¬ 
terior point of the domain is moved by V from the artificial time t = 0 to t = 1. 
At each time t , a deformation T(x,y, z,t) is constructed by the motion. Moser 
and Dacorogna proved that the final deformation, at t = 1, has the prescribed 
Jacobian determinant, i.e. its Jacobian determinant J(T(x, y, z, 1) = f(x,y,z). 
Since f(x, y , z) is positive, T(x, y , z, 1) has the desirable property of having pos¬ 
itive Jacobian determinant everywhere in the domain. 

3.2 Adaptive Non-folding Mesh Generation and Adapta¬ 
tion on Moving Domains 

We now present our deformation method on moving domains. To generate a 
mesh on a domain D 2 , we deform an initial mesh on 12i by a discrete approxima¬ 
tion of a smooth and invertible deformation <fi from 12 1 to U 2 . In fact, a family 
of adaptive non-folding meshes can be numerically generated if a positive, time 
dependent size function /(%, t) > 0 is specified, where x = ( x ,y) in 2D, t > 0 is 
a parameter (it could be the real time in unsteady field simulation). 

In this subsection, we describe the deformation method for generation of a 
family of deformations according to a prescribed (time dependent) size function 
/(%, t) > 0. We will explain how to construct higher order meshes below. 

The problem to solve is the following: Given a size function f(\, t) > 0, find 
deformations <f>(x,y,t) such that the Jacobian determinant of the deformations 
are equal to the given size function /(x, t).This problem is solved in the following 
two steps. 

Step 1: Find a vector field U by solving the divergence-curl system: 

di = 


cur4 U(<f>, t) = 

on 12i, with the Dirichlet condition on moving portion of the boundary, and the 
Neumann condition on fixed boundary. For later use we have introduced g = y- 


d 


dt 

d fl , 
at 9W>t) 
0 
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Step 2: Define the velocity vector field by V = fU and solve the ordinary 
differential equation for the transformation: 


^ 

= V{(j>,t) 

with the initial condition <j>(x,y, 0) = (a :,y). From [Cai et al. 2004] [7] we have 
the following theorem. 

Theorem 1. The deformation constructed by the above steps satisfies, Vf > 
0, det'V4>(x,t) = /(</>, t). 

Indeed, defining H = we can show ^ = 0 directly based on the 

two steps, which then implies that H is constant. For the details, we refer to 
the dissertation by Dion Fleitas. 

4 The Numerical Implementation 

We use the example in section 2 to describe the numerical implementation of 
the deformation method. 

We denote the background domain by 0(0). We will deform the initial 
Cartesian mesh on 0(0) by a suitable velocity vector from t = 0 to t = 1, and 
denote the deforming domain by 0(1). In the example, O(t) = 0(0). Suppose 
we want Pi on 0(0) to be moved to Qi on 0(1). For instance, in the example, 
we want node Pi = (1,4) to be moved to Qi = #1; P 2 = (1,2) to Q 2 = #3; 
..., Pig = (1,11) to Qi 8 = #18, etc. A pseudocode for selecting the nodes Pi as 
a general case is provided in the Appendix. 

In order to define a correct velocity vector held V, we first find a vector 
held U((f), t) by solving the divergence-curl system on the deforming domain 
0(f) at each time t = where n is the number of time steps for solving the 
ordinary differetial equation (ODE) from t = 0 to t = 1, k = 1,2, ...,n — 1. In 
the example, we take n = 10. The following equations are solved by the least 
squares finite element method (as described in [Cai et al. 2004] [7j and [Liao et 
al. 2009] [8]) at each hxed time t = ^ : 

div ‘ w - f) = hm 

= IsOM) 

curl ^ [/(</>, t) = 0 

with the Dirichlet condition U = = 1,2, 3,... ,18 in the example, and 

n-j 

U = 0 at all other boundary nodes of the background domain, where PiQi is a 
vector from Pi to Qi. Solving the ODE: 
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After determination of JJ at t = k = 1, 2 ,n — 1, we define the velocity 
vector field by V{(j), t) = f(cj>, t)U(cj), t) and solve the ODE for the transformation 
(j){ x, t) from i = 0 to f = 1, inn time steps: 


d(j)(x,t) 

dt 




with the condition <f>(x, 0) = x. 

The deformed mesh on the background domain is now divided into two 
parts: the exterior part is disregarded as shown in Figure [2a| the interior part is 
reconnected to be a triangular mesh in a straight forward manner as shown in 
Figure 2b If an exterior mesh is desired, then we disregard the interior mesh. 


A Selection of Background Mesh Nodes Pf s 

We can lay a background mesh of step size h over a contour such that each cell 
of the mesh contains only one contour point (provided the mesh won’t be too 
dense). Then we need only to find the shortest distance between the contour 
point of interest and the four points that make up the cell containing that con¬ 
tour point to avoid searching for the shortest distance between the contour point 
and all mesh points. To do this we search for the minimium distance from one 
contour point to another, say, d c , and then let h = °'^ c be the step size of the 
mesh. Let y k = {t k , s k ) be the fcth point on the contour and x- h j a mesh point. 
Let B Xi . and B Vk be elements of Boolean matrices determining if a mesh point, 
Xi j, has been moved to a contour point or if a contour point, y k , has had a mesh 
point moved to it, respectively. Then the following pseudo-code selects a set of 
background mesh points that will be moved to the corresponding contour points. 

for (k = 1 to n) 

+ tp • 

1 ~ h ’ 

* _ h ’ 

* = L* J; 
j = IaJ ; 

di = distance(xij 1 yk)', 
d 2 = distance(x i+ ij, yk)', 
d 3 = distance(x itj+ i,y k ); 
di = distance(x i+1 j +1 ,y k ); 
dmin = min{di,d 2 ,d 3 ,di)\ 

If {dmin —— di) 

if B Xi;j == false 

%i,j = V k 5 
B Xij = true ; 

B Vk = true-, 

else if ( B Xi . == true ) 
dmin = min{d 2 ,d 3 ,di)-, 
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if ( drain -6^2) 

if B Xi+1 J == false 

— Vki 

B Xi+lti = true; 

B yk = true-, 

else if (B Xi+1 . == true) 
dmin = min(d 3 , d 4 ); 
if {djnin —— d 3 ) 

if B Xi j+1 == false 

Uki 

B XiJ+1 = true-, 

B Vk = true; 

else if ( B Xi , +1 == true) 
dmin — min(d4); 
if (dmin ^4) 

if B x . +lj+1 == false 
%i+l,j-\-l = Vk\ 

B Xi+1 , j+1 = true; 

B Vk = true; 

end 
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